Earthquake & fault visualization with Leaflet for R

Jack Colpitt

The leaflet map below will be leveraging geopspatial libraries in the R ecosystem to display earthquakes within the last 30 days. Faults will also be displayed to give further insights to the clustering of seismic events.

Import libaries

library(sf) # simple features for geometries          
library(rgdal) # OGR GeoJSON driver for importing data      
library(leaflet) # mapping visualization

Data Sources / Create Spatial Dataframes

The earthquake data is being pulled in through the real-time feed provided by the USGS https://www.usgs.gov/programs/earthquake-hazards/data.

Fault polylines are provided by esri in the ArcGIS Online (AGOL) living atlas

AGOL - https://pennstate.maps.arcgis.com/home/item.html?id=37a384d4c1ef4f56a33a40f291a634e9

REST Page - https://services.arcgis.com/jIL9msH9OI208GCb/ArcGIS/rest/services/Active_Faults/FeatureServer/0

url <- "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_month.geojson" # read geojson url into a variable
earthquakes <- readOGR(url, verbose=FALSE) # use the rgdal library to read the url into a SpatialPointsDataFrame
eqsf <- st_as_sf(earthquakes) # convert sp data into a simple features dataframe

# use the query tool in the REST page (set where clause to '1=1' to get all results) to generate a geojson request url
faults_url <- "https://services.arcgis.com/jIL9msH9OI208GCb/ArcGIS/rest/services/Active_Faults/FeatureServer/0/query?where=1%3D1&objectIds=&time=&geometry=&geometryType=esriGeometryEnvelope&inSR=&spatialRel=esriSpatialRelIntersects&resultType=none&distance=0.0&units=esriSRUnit_Meter&relationParam=&returnGeodetic=false&outFields=*&returnGeometry=true&featureEncoding=esriDefault&multipatchOption=xyFootprint&maxAllowableOffset=&geometryPrecision=&outSR=&defaultSR=&datumTransformation=&applyVCSProjection=false&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnExtentOnly=false&returnQueryGeometry=false&returnDistinctValues=false&cacheHint=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&having=&resultOffset=&resultRecordCount=&returnZ=false&returnM=false&returnExceededLimitFeatures=true&quantizationParameters=&sqlFormat=none&f=pgeojson&token="

faults <-readOGR(faults_url, verbose=FALSE) # use the rgdal library to read the url into a SpatialPointsDataFrame
faults_sf <- st_as_sf(faults) # convert sp data into a simple features dataframe

Render leaflet map

The code below strings together multiple functions to render a complete web map with basemap options, class break symbology, pop-ups, layer toggling, and a legend. Lets walk through the steps…

- Set the R chunk options to ```{r, fig.width=9.5,fig.height=6} to display the width of the page.

- Use the colorBin function to create the class break symbology.

- Set the view.

- Add basemap tiles from CARTO and esri.

- Add Data.

- Configure pop-ups for earthquakes.

- Add legend and layer control.

# create class breaks with "Spectral" based on earthquake magnitude
pal <- colorBin(
  palette = "Spectral",
  domain = eqsf$mag, # use magnitude variable
  reverse = TRUE, # reverse the color direction
  bins = 5 # 5 breaks
)

leaflet() %>%
  setView(-117.841293, 46.195042, 3) %>% # set view to greater North America
  addProviderTiles(providers$CartoDB, group = "Grayscale") %>% # add CARTO tiles
  addProviderTiles(providers$Esri.WorldTerrain, group = "Terrain") %>% # add esri tiles
  addPolylines(data = faults_sf,
               popup = paste0( # create custom pop up
                 "<strong>Name:</strong> ", faults_sf$name,
                 "<br><strong>Slip Type:</strong> ", faults_sf$slip_type
                 ),
               group = "vectorData") %>% # add fault lines
  addCircleMarkers(data = eqsf, # create circle markers from earthquake data
                   fillColor = ~pal(mag), 
                   radius = ~eqsf$mag * 2,
                   stroke = FALSE,
                   color = "Spectral",
                   fillOpacity = 0.6,
                   popup = paste0( # create custum pop-ups
                     "<strong>Title:</strong> ", eqsf$title,
                     "<br><strong>Magnitude:</strong> ", eqsf$mag,
                     "<br><strong>Intensity:</strong> ", eqsf$mmi,
                     "<br><strong>Significanceq:</strong> ", eqsf$sig
                   ),
                   group = "vectorData") %>% 
  # add legend to the bottom left
  addLegend(pal = pal, values = eqsf$mag, position = "bottomleft", title = "Magnitude") %>%
  # create a layer toggle for the basemaps and vector data
  addLayersControl(overlayGroups = c("vectorData"), baseGroups = c("Terrain", "Grayscale"))

The options above from leaflet provide an interactive map with informational pop-ups for the recent seismic events and faults around the globe.